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ABSTRACT 

We present N-body simulations of resonant planets with inclined orbits that show chaotically evolv¬ 
ing eccentricities and inclinations that can persist for at least 10 Gyr. A wide range of behavior is 
possible, from fast, low amplitude variations to systems in which eccentricities reach 0.9999 and incli¬ 
nations 179.9°. While the orbital elements evolve chaotically, at least one resonant argument always 
librates. We show that the HD 73526, HD 45364 and HD 60532 systems may be in chaotically-evolving 
resonances. Chaotic evolution is apparent in the 2:1, 3:1 and 3:2 resonances, and for planetary masses 
from lunar- to Jupiter-mass. In some cases, orbital disruption occurs after several Gyr, implying the 
mechanism is not rigorously stable, just long-lived relative to the main sequence lifetimes of solar-type 
stars. Planet-planet scattering appears to yield planets in inclined resonances that evolve chaotically 
in about 0.5% of cases. These results suggest that 1) approximate methods for identifying unstable 
orbital architectures may have limited applicability, 2) the observed close-in exoplanets may be pro¬ 
duced during the high eccentricity phases induced by inclined resonances, 3) those exoplanets’ orbital 
planes may be misaligned with the host star’s spin axis, 4) systems with resonances may be systemat¬ 
ically younger than those without, 5) the distribution of period ratios of adjacent planets detected via 
transit may be skewed due to inclined resonances, and 6) potentially habitable planets in resonances 
may have dramatically different climatic evolution than the Earth. The GAIA spacecraft is capable 
of discovering giant planets in these types of orbits. 


1. INTRODUCTION 

Exoplanetary systems with multiple planets show a 
wide variety of orbital behavior , such as large amplitude 


Barnes & Greenberg 

2006 

), mean motion resonances 

(MMRs) (e.g. Lee & Beale 

2002 

Raymond et al. 20081), 

and, in one case, oscillations of inclination (McArthur 


et al. 2010 

planets in MlVIRs with mutual inclinations and find that 
these systems can evolve with large, chaotic fluctuations 
of eccentricity and inclination for at least 10 Gyr, but 
the MMR is maintained throughout. 

In general, orbital interactions are broken into two 
main categories: secular and resonant. The former treats 
the orbital evolution as though the planets’ masses have 
been spread out in the orbit, i.e. the mass distribution 
averaged over timescales much longer than the orbital 
periods. If the orbits are non-circular or non-coplanar, 
the gravitational forces change the orbital elements peri¬ 
odically. 

MMRs occur when two or more orbital frequencies are 
close to integer multiples of each other. The planets pe¬ 
riodically reach the same relative positions, introducing 
repetitive perturbations that can dominate over secular 
effects. The combination driv es the long-term ev o lution 
of the system. For example, Rivera & Lissauer (2001|) 
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considered the long-term behavior of G1 876 (Marcy et al. 


2001) and their integrations show hints of chaotic evoiu- 
tiorTon 100 Myr timescales (see their Fig. 4b). In this 
study we expand the research on exoplanets in MMRs to 
include significant mutual inclinations. 

While most bodies in our Solar System are on low-e, 
low-* orbits, the exoplanets do not share these trai ts. Ec¬ 


cent ricities span values from 0 to > 0.9 ( e.g. Butler et al. 


20061, and at least the two outer planets of u And hav e 


a large mutual inclination of 30° (McArthur et al.[2010). 


Numerous exoplanetary pairs in MiVlR are known, in- 


in 2:1 (Marcy et al.[2001 |Mayor et al.|2004 

Tinney et al. 

2006 

Vogt et al. 2005), HD 45364 in 3:2 ( 

Correia et al. 

2009 

and HD 60532 in 3:1 ( 

Desort et al. 2008). The or- 


with HST (Benedict et al. 2002]), and it was reported that 


this observation was compatible with a mutual inclina- 


tion between orbital planes (Rivera & Lissauer 2003) 
but no formal publication demonstrated that suggestion. 
The HD 128311 planetary system lies close to the 2:1 res¬ 
onance, and recent astrometric measurements have mea¬ 
sured t he orbital plane of HD 1 28311 c, but not the other 


planet (McArthur et al. 2014). Without knowledge of 
both orbital planes, the mutual inclinations are unknown 
and use of the minimum masses and coplanar orbits may 
not reveal the true orbital behavior, so this system could 
in fact lie in the 2:1 resonance. 

Several studies have also explored MMRs with mu¬ 
tual inclinations in the context of planet formation. 


Thommes & Lissauer (2003) found that convergent mi- 
gration in a planetary dis k can excite large i nclina tions 


in exoplanetary systems. Lee & Thommes (2009) per¬ 
formed a similar experiment and found that migrat¬ 
ing planets could become temporarily captured in an 











































































2 


inclination-type MMR, in which conjunction librates 
about the midpoint between the longitudes of ascending 
node (inc lination resonances are de scribed in more detail 
in § [2|. Libert & Tsiganis (2009) explored the forma¬ 
tion of higher-order inclination resonances and also found 
that temporary capture in an inclination resonance can 
occur. They also artificially turned off migration shortly 
after capture an d found the res ulting system was stable. 


Teyssandier & Terquem (20141 considered this process 
and identified several constraints on the planetary mass 
ratio and eccentricities that permit entrance into an incli¬ 
nation resonance. All these studies find that MMRs with 
inclination can form in the protoplanetary disk, but only 
for a few specific scenarios, and even then the inclination 
resonance is likely to be fleeting. None of these studies 
considered the long-term evolution of the resonant pairs. 

Significant mutual inclinations may also be formed via 
gravitational scatterin g events th at typica lly result in the 
ejection of one planet (jMarzari fe Weidcnschilling 2002 


the fundamental plane. When a system is formed, if con¬ 
junction initially occurs near one of these stable points, 
and if there is a commensurability of mean motions such 
that the conjunction longitude varies slowly, then con¬ 
junction will tend to librate. Furthermore, because the 
orbits are farthest apart at these libration centers, res¬ 
onances can reduce the likelihood of close encounters, 
further maintaining long-term orbital stability. 

Both e-type and Atype resonances are observed in our 
Solar System, but e-type is far more prevalent, including 
the Galilean satellite system and the Neptune-Pluto pair. 
An i-type resonance is obse rved in the Saturnian satellit e 
pair of Mimas and Tethys ( Allan|1969 Greenberg|1973 ). 


Neptune and Pluto are particularly re levant to this study. 


as the pair is in the 3:2 e-resonance (Cohen & Hubbard 


19651, and later studies found that the i-resonance ar- 


Chatterjee et al. 

2011:). R 


yrriond et al.|2010||Barnes et al. 

wed tna 


L|). [Raymond et al. (|2C)10|) also showed tnat scattering 

a produce systems in MMRs about 5% of the time, 


guments also librate on short timescales, but the libra- 
tion drifts slowly such that the resonant arg ument actu ¬ 
ally circulates with a period of ~ 25 Myr (Williams & 
Benson|[l971 Applegate et al. 1986). Higher order secu- 


coulc 

but they did not consider the inclinations of the resul¬ 
tant systems. All these studies reveal that formation of 
MMRs with mutual inclination is possible, but probably 
rare. 

While HST has successfully characterized several 
nearby exoplanetary systems astrometrically, the GAIA 


lar resonance s are also o peratin g on Neptune and Pluto 
( Milani et al.||198~9 Kinoshita fe Nakai]|1996 ), and likely 


contribute to their configuration being formally chaotic 
(Sussman fe Wisdom] If 


No ne theless, t he pai r is sta¬ 
ble tor at least 5.5 Gyr ( Rmoshita fe Nakai|1996 ). Thus, 


Pluto and Neptune’s orbital evolution is impacted by the 
i-resonance but it is likely a small effect. 

The first step in identifying a mean motion resonance is 


space telescope could astrometrically detect hundreds 

to examine the ratio of orbital periods P or, equivalently, 

of Jupiter-sized planets (see e.g. 

Lattanzi et al. 2000 

the mean motions n = 2ir/P. If the ratio is close to a 

Sozzetti et al. 2001; Casertano et a 

. 200§[ Sozzetti et al. 

ratio of two integers, i.e. ni /712 ~ ji/j -2 where j\ and 

2014), revealing how common are systems in an MMR 

with significant mutual inclinations. Against this back¬ 
drop, we explore the dynamics of planets with orbital 

j '2 are integers, then the system may be in resonance. 
An MMR requires a periodic force to applied near the 
same longitude relative to an apse or a node, which pre- 


period commensurabilities, significant eccentricities and 
mutual inclinations. 

This paper is organized as follows. In § [2] we describe 
the physics of resonance and our numericalmethods. In 
§ [3j we present results of hypo thetical systems, including 


planets in the hab itable zone (Kasting et al.||1993| Kop- 
|parapu et al.|2013 ), as well as giant and dwarf planets m 


cess with time. To account for this evolution, celestial 
mechanicians have introduced the resonant argument, a 
combination of angles that account for both the mean 
motions and the orbital orientations. 

For e-type resonances, the resonant arguments are 


#1,2 = J 1 A 2 — J 2 A 1 — J3?x7l,2, 


( 1 ) 


the 2:1, 3:1 and 3:2 commensurabilities. In §[f]we show 
that planet-planet scattering can produce systems in the 
2:1 MMR and with significant mutual inclinations. In 
§ [5] we analyze several known exoplanetary systems and 
find several that could be evolving chaotically. In § [6] we 
discuss the results and then conclude in § [7| 


2. METHODS 
2.1. Resonant Dynamics 

MMRs are well-studied, and can be explained intu¬ 
itively forjoiwbutjiomzero, va!ues_of_e_jiidJnclination i 
(|Peaie||1976| |Greenberg||1977| | Murray & Dermott||1999[) . 
Stable resonances can be divided into two types: eccen¬ 
tricity and inclination. The difference lies in the loca¬ 
tions of the stable longitudes of conjunction, sometimes 
called the libration center. In an eccentricity (e-type) 
resonance, the stable longitudes are located at the lon¬ 
gitude of periastron of the inner planet (wi), and the 
apoastron of the outer planet («2 = ^2 + 7r). For in¬ 
clination (i-type) resonances, the stable longitudes lie 
halfway between the longitudes of ascending node, Q, 
of each planet (fli ,2 ± 7r/2) when the reference plane is 


where the j terms are integers that sum to 0 , and the sub¬ 
scripts 1 and 2 to vj correspond to the inner and outer 
planet, respectively. Should any 8 librate with time, or 
even circulate slowly, then resonant dynamics are impor¬ 
tant. Conjunction (ji \2 — J 2 A 1 ) lies close to a particular 
longitude relative to the apsides. Usually 8 will librate 
about either 0 or 7 r, but other stable values are possible 
and are referred to as asymmetric resonances. 

For i-type resonances, the dominant resonant argu¬ 
ment is 

<t> = J 1 A 2 — J 2 A 1 — j3^l — J4^2) (2) 

if i = 0 corresponds to the invariable, or fundamen¬ 
tal, plane, i.e. the plane perpendicular to the total 
angular momentum of the system. With this choice, 
fix = U 2 ± 7 r. Inclination resonances are fundamen¬ 
tally different from e-type in that first order resonances 
are technically not possible, i.e. j\ — J 2 > 1- The sta¬ 
ble conjunction longitudes correspond to the two mid¬ 
node longitudes. Note that for circular orbits, there is 
no substantive difference between the two stable longi¬ 
tudes: The geom etry at ± 90° from the mutual node are 
identical. See Greenberg (1977) for more discussion on 
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the physics of the i-resonance. 

It is often convenient to think of resonances in terms 
of a pseudo-potential. The conjunction longitude, A c is 
accelerated toward the stable points, and oscillates sinu¬ 
soidally about it. The acceleration of A c in a 2:1 (4:2) 
MMR can be written as 

A c = Ciei sin(A c —n 7 i)-|-C 2 e 2 sin(ra 72 —A c )+C 3 ^ sin 2(A C —IR), 

( 3 ) 

where Ci, C 2 and C 3 are con stants that depe nd on the 
masses and semi-major axes (Greenb erg|1 977), and A c = 
2 A 2 — Ai. Eq. <[3j) is analogous to that of a compound 
pendulum, as long as the e’s and V s are approximately 
constant. The pseudo-potential contains minima at w 1 , 
02, and ili ± 7r/2, but each has different depths which 
vary as the orbits evolve. Moreover, as e and/or i become 
large, this simple picture breaks down and more minima 
will likely appear. Thus, we should anticipate the motion 
to become complicated, an expectation that is borne out 
in §§H] [5] 

2 .2. Numerical Methods 

The classical analytic theory summarized above be¬ 
comes less accurate as the values of eccentricity and in¬ 
clination increase. To lowest order, the evolution of e 
and i are decoupled, but as either or both become large, 
pathways for the exchange of angular momentum open 
up, and the m otion can become very complicated ( e.g. 
Barnes et al. 2011). Thus, we rely on IV-body numer- 
ical methods to analyze resonances with mutual incli¬ 
nations, while appealing to published analytic theory to 
help interpret the outcomes. In particular we us e the 
well-test ed and r eli able Mercury (Chambers 1999) and 
HNBody (Rauch & Hamilton 2002) codes We do not 
include general relativistic corrections, which are mostly 
negligible for the planetary systems we consider here. We 
used both mixed variable symplectic and Bulirsch-Stoer 
methods to integrate our systems. While the former 
can integrate our hypothetical systems to 10 Gyr within 
about 10 days on a modern workstation, the Bulirsch- 
Stoer method, which is more accurate, requires about 80- 
90 days with HNBody; hence we used it sparingly. Regard¬ 
less of our choice of software and/or integration scheme, 
identical initial conditions produced qualitatively similar 
results. All integrations presented here conserved energy 
to better than 1 part in 10 , usually better by many or¬ 
ders of magnitude. Unless stated otherwise, all systems 
maintained an MMR for 10 Gyr and met our energy con¬ 
servation requirements. The reference plane for i and H 
for all cases is the invariable plane so that the inclina¬ 
tions and longitudes of ascending node are more physi¬ 
cally meaningful Here we ignore planetary and stellar 
spins: Planets and stars are point masses. 

We consider several broad categories of exoplanetary 
systems. First we model systems in the 2:1 resonance. 

In Set #1, one planet is always 1 Earth-mass (1 M 0 ) 
with a semi-major axis a of 1 AU, and the pri mary is a 
solar-mass s tar, i.e. it is in the habi table zone (|Kasting 
et al. 1993 Kopparapu et al. 20131. The other planet 
may have a mass between 3 and 40 Mg. In Set #2, we 


6 We have made our source code to rotate any astrocen- 
tic orbital elements into the invariable plane, as well as gen¬ 
erate input files for Mercury and HNBody, publicly available at 
https: //github.com/RoryBarnes/InvPlane. 


explore the case of two giant planets orbiting a 0.75 M 0 
star in 2 and 4 year orbits. Set #3 considers two dwarf 
planets orbiting a solar-mass star. In Sets #4 and #5, we 
simulated ~Earth-mass planets in the 3:1 and 3:2 MMR, 
respectively. For all these simulations, the period ratios 
are at exact resonance, but the other orbital elements are 
chosen randomly and uniformly over a wide range of val¬ 
ues. For each set, we integrated 100 systems for 10 Myr 
as an initial survey, and for systems that appeared inter¬ 
esting, i.e. showed chaotic evolution, we continued them 
to 10 Gyr using 0.01 year timesteps. Table [T] shows the 
ranges of initial conditions we used for these sets of hy¬ 
pothetical systems. 

In § 5 we examine known systems in or near the 2:1 
(HD 73526 and HD 128311), 3:2 (HD 60532) and 3:1 
(HD45364) MMRs. Each of these systems was discov¬ 
ered via radial velocity data, and hence suffer from the 
mass-inclination degeneracy. However, HD 12 8311 c has 
also been d etected astrometrically with HST ( McArthur] 
et al.|20141, hence its mass and full orbit are known. The 
best fits and uncertainties in parentheses for these known 
systems are listed in Table [2] The position of each planet 
in its orbit, the “phase,” is crucial information and differ¬ 
ent authors present the information in different formats 
and so we present the parameter used in the most recent 
paper in Table 2. T p is the time of periastron passage, A 
is the mean longitude and /x is the mean anomaly. 

For each of the known systems, we vary each orbital 
element uniformly within its published uncertainties and 
simulate the orbital evolution with HNBody. Our goal is 
not to calculate a probability that each system is in an 
i-type resonance, but rather to determine if it is possible 
at all. We simulate 100 versions of the these systems 
and allow i and H to take any value, and m is adjusted 
accordingly. If chaotic motion in the MMR is apparent, 
we integrate it for a further 10 Gyr. 


3. HYPOTHETICAL SYSTEMS 

In this section we present the results of the simulations 
described in the previous section. We separate the results 
by resonance: 2:1, then 3:1, and finally 3:2. In all cases 
we find evidence of chaotic orbital evolution, often with 
very large amplitudes of eccentricity and inclination. 


3.1. The 2:1 Resonance 

In this section we examine example cases in Sets #1 3, 
i.e. in the 2:1 MMR. The e-resonance arguments are 

6 * 1,2 = 2 A 2 — Ai — K7i,2) (4) 

and the i-resonance argument is 

cf> = 4 A 2 — 2Ai — Hi — H 2 . (5) 

3.1.1. Earth with an Exterior Companion (Set #1) 

Set #1 consists of an Earth-mass planet with a 1 year 
period and a larger exterior planet with an orbital period 
of 2 years. In Fig.[l]we show the evolution of the resonant 
arguments, e ana! for the first example, System A in 
Table [ 3 ] for the first 10 5 years. 

The top two panels show that the resonance arguments 
switch between libration and circulation. Initially 0\ 
(black dots) librates about 0, the classic libration cen¬ 
ter for e-type resonances, while 62 (red dots) librates on 
short timescales, but circulates on long timescales. Also 
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TABLE 1 

Initial Conditions for Hypothetical Systems 


Set 

MMR 

M* (M 0 ) 

mi 

ra 2 

oi (AU) 

a 2 (AU) 

e 

<(°) 

n(°) 

" (°) 

A(°) 

i 

2:1 

1 

1 M@ 

3.3-36.6 M® 

1 

1.5874 

0-0.5 

0-30 

0-360 

0-360 

0-360 

2 

2:1 

0.75 

0.314 M Jup 

0.1 1.15 M Jup 

1.44219 

2.2893 

0-0.5 

0-30 

0-360 

0-360 

0-360 

3 

2:1 

1 

0.0775 M Moon 

2.84 M Moon 

1 

1.5874 

0-0.5 

0-30 

0-360 

0-360 

0-360 

4 

3:1 

1 

1 M® 

3.3-36.6 M® 

1 

2.0801 

0-0.5 

0-30 

0-360 

0-360 

0-360 

5 

3:2 

1 

1 M® 

3.3-36.6 M® 

1 

1.3104 

0-0.5 

0-30 

0-360 

0-360 

0-360 


TABLE 2 

Best Fits and Uncertainties for Selected Known Exoplanet Systems 

System 

M, (M 0 ) 

Planet 

m (M Jup ) 

P(d) 

e 

" (°) 

Phase 

HD 128311 

0.828 

b 

1.769 (0.023) 

453.019 (0.404) 

0.303 (0.011) 

57.864 (3.258) 

2400453.019 (4.472)“ 



c 

3.125 (0.069) 

921.538 (1.15) 

0.159 (0.006) 

15.445 (6.87) 

2400921.538 (18.01)“ 

HD 73526 

1.08 

b 

2.8 (0.2) 

188.3 (0.9) 

0.19 (0.05) 

203 (9) 

86 (13) 6 



c 

2.5 (0.3) 

377.8 (2.4) 

0.14 (0.09) 

13 (76) 

82 (27 ) b 

HD 60532 

1.44 

b 

1.03 (0.05) 

201.3 (0.6) 

0.28 (0.03) 

351.9 (4.9) 

2453987 (2)“ 



c 

2.46 (0.09) 

604 (9) 

0.02 (0.02) 

151 (92) 

2453723 (158)“ 

HD 45364 

0.82 

b 

0.1872 (0) 

226.93 (0.37) 

0.1684 (0.019) 

162.58 (6.34) 

105.76 (1.41) c 



c 

0.6579 (0) 

342.85 (0.28) 

0.0974 (0.012) 

7.41 (3.4) 

269.52 (0.58)“ 


Tperi (JD) 

b M (°) 
c A (°) 


TABLE 3 

Initial Conditions for Selected Set #1 Systems 


System 

Body 

m (M®) 

a (AU) 

e 

i (°) 

« C) 

u C) 

M(°) 

A 

1 

1 

1 

0.0866 

10.322 

223.75 

340.2 

268.33 


2 

10.07 

1.5874 

0.1883 

0.82 

43.747 

74.21 

296.5 

B 

1 

1 

1 

0.0251 

1.61 

277.19 

353.67 

261.33 


2 

26.91 

1.5874 

0.0531 

0.0475 

97.19 

219.07 

359.38 

C 

1 

1 

1 

0.2373 

2.918 

217.13 

330.25 

102.45 


2 

4.39 

1.5874 

0.4225 

0.565 

37.13 

246.36 

112.88 

D b 

1 

1 

1 

0.00296 

19.83 

278.01 

332.27 

268.52 


2 

35.62 

1.5874 

0.266 

0.449 

98.01 

69.73 

343.09 

E 

1 

1 

1 

0.2952 

38.51 

247.92 

40.44 

225.93 


2 

14.82 

1.5874 

0.0961 

1.832 

67.92 

253.92 

318.12 

F 

1 

1 

1 

0.1 

12.04 

270.26 

140.77 

0 


2 

10 

1.5874 

0.15 

0.9551 

90.26 

23.77 

10 

G 

1 

1 

1 

0.15 

23.526 

232.96 

31 

0 


2 

10 

1.5874 

0.2 

1.832 

52.96 

247.9 

10 


b Stable for only 73 Myr. 


note that the e- and ^-resonance arguments appear to 
be coupled. During this initial phase, <f> librates with 
large amplitude. After about 14,000 years, the behav¬ 
ior changes dramatically and all arguments librate, but 
with sudden jumps between libration centers. Then at 
25,000 years, the motion appears to return to the ini¬ 
tial state. This switching between modes of oscillation 
demonstrates the system is chaotic, and is often an indi¬ 
cator of impending instability (e. g. Laskar|1 990). Hence 
we would naively expect a similar outcomeTor this sys¬ 
tem. The eccentricity of the inner planet shows un¬ 
usual behavior that also changes with the resonance ar¬ 
guments. The evolution is approximately periodic, but 
does not appear regidar. The inclinations do not appear 
to be strongly impacted by the changes in the resonance 
argument behavior. 

In Fig. [2] we extend the evolution of e and i by a factor 


of 10, to 10 Gyr. Here we see that the hints of chaotic 
behavior in Fig. |T] remain present, but are not indica¬ 
tive of the true scale of the chaotic motion for this sys¬ 
tem. Remarkably, the system survives for 10 Gyr despite 
the chaos. The eccentricity of the inner planet aperiodi- 
cally reaches values larger than 0.65, while its inclination 
reaches 40°. The outer planet is more massive, so its evo¬ 
lution does not have as large an amplitude as the inner, 
but it, too, evolves chaotically. While the variations in 
e and i are chaotic, there are clearly bounds to their 
permitted values. 

The libration and slow circulation of the resonant ar¬ 
guments suggest that both e and i-type resonances are 
important in this system. We further explore the evo¬ 
lution of the system in Fig. |3l which only considers the 
first 30,000 years of orbital evolution. In the left panel we 
compare the conjunction longitude (black dots) with the 
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0 2-10 4 4*10 4 6*10 4 8*10 4 MO 5 

Time (yr) 


Fig. 1.— The first 10 5 years if evolution of System A. Black corresponds to the inner Earth-mass planet initially at 1 AU, red to a larger 
planet in the outer 2:1 resonance. Variations of the inclination resonance argument (top), eccentricity arguments (top middle) with black 
dots corresponding to the zu\ argument and red to vjz, eccentricities (bottom middle) and inclinations (bottom). 
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Fig. 2.— The evolution of e and i for System A over 10 Gyr. Black points correspond to the inner planet and red to the outer. 
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various stable longitude predicted from classic resonant 
theory, w\ (purple curve), a 2 = ru 2 +7r (orange curve), 
flf = f2i + 7r/2 (blue curve), and — n/2 

(green curve). Conjunction librates, but, surprisingly, 
not always about the expected stable longitudes. For 
most of the time, A c ~ wx, but it also librates about 
other locations that are not associated with any classic 
libration center. Note that conjunction avoids a 2 . 

In the right panel, we plot the two mean motions and 
see that they librate about the resonant frequencies, but 
with varying amplitudes. Careful inspection of the two 
panels shows that the different mean motion amplitudes 
correspond to the different librations seen in the left 
panel. There appear to be 5 different resonant oscilla¬ 
tions over this cycle, which approximately repeats for at 
least the first I Myr. 

The long-term behavior in Fig.[2]shows mode-switching 
can occur on longer timescales as well. For example, from 
6.5 - 6.8 Gyr, e and i are confined to narrow regions, but 
then return to the large amplitude oscillations. Fig. [4] 
shows the evolution over 10 5 years starting at 6.75 Gyr 
within this alternative mode. In comparison to Fig. [T] 
the ^-resonance argument is circulating, as is 9 2 . How¬ 
ever, 9\ is librating about 0 indicating that the system 
is exclusively in an e-type resonance. 

After 10 Gyr the evolution is qualitatively similar to 
the initial evolution. The e- and i-resonant arguments 
are switching between libration and circulation, a quasi¬ 
stable behavior. We are unaware of any study that has 
found qualitatively similar behavior in a planetary or 
satellite system. W hile lon g-li ved chaos is evident in our 


Solar System ( e.g. Sussman & Wisdom 
plitudes of the variations are much sma 


1988), the am- 
ier. Moreover, 


the switching between different quasi-periodic states is 
also usually an indicator of instability. System A is not 
“nearly integrable,” as is often argued for stable plane¬ 
tary systems. Nonetheless, these planets remains bound 
and in resonance for the main sequence lifetime of a solar- 
type star. 

Although the orbital behavior of System A is surpris¬ 
ing, this case is not anecdotal. Table [3 contains 6 exam¬ 
ples from Set #1 that show chaotic evolution, yet remain 
in resonance for 10 Gyr. In Figs. [5]j7 we show 3 more 
cases as representative examples ofthe 2:1 MMR. These 
simulations demonstrate that chaotic inclined MMRs can 
produce behavior that spans a range from high frequency, 
low amplitude oscillations to low frequency and high am¬ 
plitude oscillations to cases in which all available phase 
space is sampled. 

System B is an example of high frequency, low ampli¬ 
tude evolution. Fig. [5] shows the first 25,000 years of its 
evolution. In this case the resonant arguments oscillate 
with a period of a few hundred years and switch states at 
^18,000 years. The eccentricities and inclinations of the 
Earth-like planet oscillate with an amplitude of 0.1 and 
0.25°, respectively, on this timescale. Note that System 
B begins with a mutual inclination of just 1.65°, reveal¬ 
ing that relatively small mutual inclinations can lead to 
chaotic orbital evolution. Over 10 Gyr, e\ remains below 
0.14 and e 2 0.06, while i\ remains below 18° and i 2 below 
1 °. 

System C is similar to System A in its evolution, as 
shown in Fig. [7] Over 10 Gyr e\ varies from 0 to 0.9 and 


it varies from 0 to 60°. From ~ 3 Gyr to ~ 6 Gyr (not 
shown) the system enters a different mode in which the 
eccentricities and inclinations are confined to narrower 
regions. 

System D only survives in the resonance for 73 Myr, 
but we it include here to illustrate the extreme eccen¬ 
tricity evolution that is possible in inclined MMRs. At 
87,170 years, e\ reaches a value of 0.99998, implying a 
periastron distance that would place it inside the core of 
its solar-mass primary. Clearly, the point-mass approxi¬ 
mation for the orbital dynamics has broken down for this 
system. In particular, we expect that at when e\ ~ 1 
that strong tidal effects should dramatically alter the or¬ 
bits. We return to this point in § [6] 

Rather than show similar plots for Systems E-G, we 
only comment on their behavior. System E shows a cir¬ 
culating <f> for the first 10 5 years, while 9 1 librates with 
large amplitude about 0, and 0 2 drifts. The inner planet’s 
e and i vary from 0 to 0.75 and 55°, respectively. After 
~ 6 Gyr, the system slowly moves into a different state 
with ix varying between 20° and 50°, and 9\ aperiodi- 
cally circulating. System F has resonant arguments that 
switch between libration and circulation as in Fig. [l] and 
eccentricities that remain below 0.3 and inclinations be¬ 
low 20°. System G is similar to System A. 

These examples are only illustrative and do not rep¬ 
resent the full range of possibilities. The 7 cases listed 
in Table [3] met our stability criteria (see § [ 2 ]) , and we 
suspect many more cases also would, but computational 
constraints prevented a more thorough analysis. 


3.1.2. Two Giant Planets Orbiting a 0.75 Mq Star (Set #2) 


Next we consider Set #2: Two Saturn- to Jupiter-mass 
planets orbiting a 0.75 Mq star in 2 and 4 year orbital 
periods. Such planets induce a larger astrometric signal 
in the host star, and are large enough that radial velocity 
observations could break the 180° ambiguity in Q implicit 
in astrometric measurements of exoplanets. In Table [4] 
we present 7 systems which survived for ~ 10 Gyr and 
conserved energy adequately. 

In Fig. [8] we show the evolution of System N on two 
timescales. The left panels show the variation in the reso¬ 
nant arguments and e and i for 10 5 years in the same for¬ 
mat as Fig. [T] As before the resonant arguments switch 
between libration and circulation leading to chaotic evo¬ 
lution of e and i. In the right panels we show the evolu¬ 
tion of e and i for 10 Gyr. e\ aperiodically reaches values 
of ~ 0.85, and i\ reaches 50°. Note as well at 2.4 and 
6.0 Gyr the systems enters qualitatively different states 
for about 100 Myr. 

In Fig. [9] the orbit of the host star about the system’s 
barycenter is shown over 7 years if viewed face-on, i. e. the 
invariable plane is parallel to the sky plane. In this ex¬ 
ample we assume the system is located at 25 pc, and the 
combined system induces a ~ 100 fi as astrometric signal, 
which is 4-5 times larger than the expected GAIA uncer¬ 
tainties for stars with G-band magnitudes <13, repre¬ 


sented by t he line labeled “GAIA Uncertainty” (Sozzetti] 


et al. 2014). This system would be relatively easy to 
characterize with radial velocity data as well, and hence 
could be fully characterized in the next 10 years. Of 
course, the details of actually modeling resonant systems 
are non-trivial (e.g. Marcy et ahpOOl I, but the discov- 
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Fig. 3.— Initial evolution of System A. Left: Evolution of the conjunction longitude \ c (black dots) in relation to the 4 classic stable 
longitudes, w\ (purple curve), 0.2 = "2 +rr (dashed orange), 1 l( = fli + 7r/2 (blue), and = 11 \ — 7r/2 (green). Right: Evolution 
of the mean motions of the inner planet (top) and outer planet (bottom). 


TABLE 4 

Initial Conditions for Selected Set #2 Systems 


System Body m ( Mj up ) a (AU) e i (°) f2 (°) u (°) p (°) 


H 

1 

0.314 

1.4422 

0.0276 

2.504 

247.16 

94.57 

7.74 


2 

0.573 

2.2893 

0.2021 

1.118 

67.12 

126.15 

313.87 

I e 

1 

0.314 

1.4422 

0.1155 

3.302 

39.32 

273.88 

195.14 


2 

0.303 

2.2893 

0.4692 

3.076 

219.29 

247.62 

48.35 

J 

1 

0.314 

1.4422 

0.2374 

7.089 

24.56 

256 

77.1 


2 

0.234 

2.2893 

0.3988 

8.019 

204.54 

185.86 

127.68 

K 

1 

0.314 

1.4422 

0.0364 

11.267 

126.11 

284.5 

41.42 


2 

0.328 

2.2893 

0.223 

8.761 

306.18 

273.49 

160.47 

L 

1 

0.314 

1.4422 

0.0263 

19.459 

176.25 

202.97 

229.91 


2 

0.736 

2.2893 

0.3696 

6.989 

356.23 

153.08 

261.74 

M 

1 

0.314 

1.4422 

0.4596 

3.282 

106.68 

44.59 

350 


2 

0.399 

2.2893 

0.3402 

1.936 

286.57 

131.7 

73.93 

N 

1 

0.314 

1.4422 

0.0561 

21.623 

247.55 

325.64 

274.59 


2 

0.902 

2.2893 

0.4466 

6.53 

67.64 

60.26 

326.55 


e Stable for only 9.761 Gyr. 


ery of such a chaotic system would mark an important 
milestone in exoplanet science. 

In Fig. [lO] we show the evolution of System I in the 
same format as Fig. [8] This system appears qualitatively 
similar to System NTbut with lower amplitudes in e and 
i. However, this system destabilizes at 9.761 Gyr, as seen 
on the right side of the right panels. We found several 
other systems in which an ejection occurred after 1 Gyr. 
The chaotic resonant behavior shown in this study is not 
necessarily stable on arbitrarily long timescales. Like our 


own Solar System, these hypot hetic al s_ystems are just 
relatively long-lived (see, e.g. Lecar et ahpOOl). 


3.1.3. Dwarf Planets in the 2:1 MMR (Set #3) 

In this section we consider Set #3, dwarf planets in 
the 2:1 MMR. While these planets are not detectable in 
the near future, they display remarkable dynamics and 
illustrate the extreme chaos that the 2:1 MMR with in¬ 
clinations can produce. In Table [5| we pre sent 4 cases 
that are stable for 10 Gyr, and in Fig. 11 we show the 
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Fig. 4.— Same as Fig.|^ but at 6.75 Gyr. 
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Fig. 5.— The first 25 kyr of evolution of System B in the same format as Fig. ^ 








11 



0 2-10 5 4*10 5 6*10 5 8*10 5 MO 6 

Time (yr) 


Fig. 6.— The first 1 Myr of evolution of System C in the same format as Fig. ^ 
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Fig. 7.— The first 10 5 years of evolution of System D in the same format as Fig.JT] This system destabilized after 73 Myr. 
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Orbital evolution of System N. The left panels are in the same format as Fig. [T] the right is in the same format as Fig. [2] 


Fig. 8.— 

evolution of System P over 10 Gyr. Note that the e- 
resonance arguments sometimes librate for long periods 
of time, such as near 2.5 Gyr. The i argument does not 
appear to librate in this visualization, but higher reso¬ 
lution plots over shorter periods show similar behavior 
as above. Systems O and Q are similar to P, but the 
inclinations and eccentricities remain lower. 

These systems show a diversity of behavior from small- 
scale chaos (System Q), to dramatic chaos in which ei 
and ii sample all available phase space (System P). The 
resonant arguments, particularly the eccentricity argu¬ 
ments, switch between circulation and libration. e\ in 
Systems O (not shown), P and R (not shown) aperiodi- 
cally reach values in excess of 0.99, and hence they should 
tidally circularize prior to 10 Gyr. 

The resonant arguments in these systems behave dif¬ 
ferently than for the larger planets. On short timescales 
(not shown) the resonant arguments switch modes as 
in previous cases, but these systems can remain in one 
mode f or l ong timescales, particularly the e-resonance. 
In Fig. [lT] note that there are intervals when the two 
e-arguments librate for more than 100 Myr. 

3.2. The 3:1 Resonance 

In this section we consider an Earth-like planet at 1 
AU from a 1 M 0 star with an exterior companion with 
an orbital period of 3 years, i.e. in the 3:1 MMR. In this 
case the e-resonance arguments are 

( 6 ) 


and the i-resonance argument is 

4> = 3 A 2 — Ai — fli — tf.2- (7) 

Table [ 6 ] lists three configurations that are stable for 
10 Gyr and show chaotic evolution. Systems T and U 
are shown in Figs. [12] and |13[ respectively. The former 
has (relatively) modest inclination and eccentricity vari¬ 
ations, while the latter samples all the phase space avail¬ 
able, with 0 <; * <5 7 r and 0 <; e <; 1. In both systems the 
resonance arguments switch between libration and circu¬ 
lation, as seen in the previous 2:1 MMR cases. System S 
(not shown) is similar to System U, but the inclinations 
only reach ~ 65° and e only 0.9. 

3.3. The 3:2 Resonance 

Next we consider the 3:2 MMR with an interior 1 M 0 
planet at 1 AU from a solar-mass star and a larger exter¬ 
nal companion at 1.3104 AU. The e-resonance arguments 
are 

0i,2 = 3A2 — 2Ai — (8) 

and the i-resonance argument is 

(j> = 6 A 2 — 4Ai — Ui — U 2 . (9) 

Tablc[7]lists two systems that are stable for 10 Gyr and 
showed chaotic evolution. In Fig. [l4]we plot the evolu¬ 
tion of the resonant arguments, e and * on two timescales 
for System W. The evolution is qualitatively similar to 
that in the other resonances. System V is quite similar, 
with e reaching 0.9 and i reaching 60° aperiodically over 
10 Gyr. 


$1,2 — 3A2 — Ai — 2 h7i j 2, 
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Fig. 9.— Orbit of the host star in Fig. 18] (System N) about the system barycenter for 7 years if the system lies at 25 pc and the fundamental 
plane is parallel to the sky plane. The nne labeled “GAIA Uncertainty” is 25 fia,s long, and represents GAIA’S typical per-measurement 
uncertainty for bright stars. 
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Fig. 10.— Orbital evolution of System I in the same format as Fig. [8] The system destabilizes after 9.761 Gyr. 


TABLE 5 

Initial Conditions for Selected Set #3 Systems 


System 

Body 

TYl (M Moon) 

a (AU) 

e 

i (°) 

n (°) 

w (°) 

A 4 (°) 

o 

i 

0.0775 

1 

0.3032 

19.835 

283.31 

290.53 

231.72 


2 

1.171 

1.5874 

0.2768 

1.013 

103.31 

50.74 

54.65 

p 

1 

0.0775 

1 

0.04305 

24.34 

57.06 

90.49 

57.09 


2 

2.817 

1.5874 

0.3192 

0.544 

237.06 

231.06 

90.68 

Q 

1 

0.0775 

1 

0.1036 

5.298 

169.03 

21.27 

348.68 


2 

0.568 

1.5784 

0.1544 

0.577 

349.03 

171.19 

158.69 

R 

1 

0.0775 

1 

0.4475 

14.98 

289.05 

94.15 

272.54 


2 

2.513 

1.5784 

0.2093 

0.332 

109.05 

352.99 

154.19 


TABLE 6 

Initial Conditions for Selected Set #4 Systems 


System 

Body 

m (M®) 

a (AU) 

e 

i (°) 

n (°) 

«(°) 

M(°) 

s 

1 

1 

1 

0.1798 

2.06 

178.04 

261.79 

66.85 


2 

6.5 

2.0801 

0.3808 

0.234 

358.4 

113.3 

286.96 

T 

1 

1 

1 

0.02514 

1.04 

246.96 

151.63 

261.33 


2 

27.05 

2.0801 

0.05311 

0.0362 

66.96 

17.04 

359.38 

U 

1 

1 

1 

0.149 

43.64 

260.51 

18.15 

220.4 


2 

22.2 

2.0801 

0.2755 

1.27 

80.51 

51.82 

6.72 
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Fig. 11.— Evolution of System P in the same format as Fig. [T] 
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Fig. 12.— Orbital evolution of System T in the same format as Fig. [8] 



TABLE 7 

Initial Conditions for Selected Set #5 Systems 


System 

Body 

m (M e ) 

a (AU) 

e 

i (°) 

n (°) 

" (°) 

M(°) 

V 

i 

1 

1 

0.1404 

27.34 

232.78 

39.79 

5.92 


2 

5.96 

1.3104 

0.3498 

4.08 

52.78 

314.14 

108.38 

w 

1 

1 

1 

0.0144 

36.05 

316.73 

13.72 

293.02 


2 

21.83 

1.3104 

0.1671 

1.37 

136.73 

307.8 

308.84 
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Fig. 13.— Orbital evolution of System U in the same format as Fig. 

4. FORMATION BY SCATTERING 

The large amplitude chaotic orbital evolution shown 
above is remarkable, but can such a system form? As 
described in § [T| several studies have examined the for¬ 
mation of inclination resonances in the context of con¬ 
vergent migration during the protoplanetary disk phase. 
Those studies found that i-resonances can form under the 
proper circumstances, but they did not consider the post¬ 
formation evolution. In this section, we show that grav¬ 
itational scattering among planets can produce systems 
in the 2:1 MMR with mutual inclinations that evolve 

chaotically for 10 Gyr. _ 

Our samp le comes from the data set used in Raymond 
et al. (2008) that found MMRs resulting from scattering 
events. The reader is referred to that paper for details, 
but briefly, systems consisting of initially 3 planets were 
allowed to interact gravitationally and in many cases 1- 
2 plan ets were ejected from the system. Raymond et al. 


(2008) examined the e-resonance arguments of systems 
in which 2 planets remained, and reported that 5% of 
systems had at least one e-resonance argument librating. 
Systems like those shown in § [3] were deemed unstable 
and thrown out. 

In light of the results of § [3] we have re-examined sys¬ 
tems near an MMR to searcn for chaotic, but long-lived 
evolution. Specifically we focus o n the “Mixed2” distri¬ 
bution of Raymond et al.| (|2008|) in which the planets’ 
masses follow a power-law distribution with exponent 
— 1.1. This distribution has recently been shown to repro¬ 
duce many observed dynamical properties of exoplanets 


(Timpe et al.||2013]). This suite of simulations consisted 
of 1000 systems, and we examined 49 that had period 
ratios with 10% of 2, of which 24 were i dentif ied as being 
in an e-resonance by Raymond et al. (20081. Of these, 
27 had mutual inclinations less than and the largest 
was 27°. We find three that appear qualitatively similar 
to those in the previous section. Th ey a re listed in Ta¬ 
ble |8j and System X is shown in Fig. [l5] System Y (not 
shown) evolved such that the eccentricities and inclina¬ 
tions remained below 0.12 and 7°, respectively. System 
Z evolved such that they remained below 0.15 and 3°. 
We also searched for systems evolving chaotically in the 
3:1 or 3:2 MMR, but did not find any. 

The amplitudes of the variations of the orbital elements 
is lower than some cases shown in § [3j but similar to 
others, e.g. System B. Given the small number of systems 
that produced chaotic resonant behavior, it remains to 
be seen if evolution that reaches e ~ 1 and * ~ n can 
be naturally produced. Nonetheless we conclude that 
planet-planet scattering can produce systems that evolve 
chaotically for 10 Gyr in an MMR. 

5. KNOWN SYSTEMS 

The previous two sections established the typical char¬ 
acteristics of planets in inclined MMRs and a viable for¬ 
mation mechanism. The next question that naturally 
arises is if any known systems might be evolving in a 
long-lived and chaotic configuration. In this section we 
examine four candidates in 3 different MMRs and with 
observed properties listed in Table [2] These planets were 
all discovered by the radial velocity technique and hence 
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TABLE 8 

Initial Conditions for Planets in a Chaotic 2:1 MMR Formed by Scattering 


System 

Body 

m (M Jup ) 

a (AU) 

e 

i(°) 

n (°) 

" (°) 

M(°) 

X 

1 

0.118 

4.83696 

0.1498 

18.845 

291.14 

83.60 

330.83 


2 

0.213 

7.71971 

0.1548 

8.15 

111.12 

195.95 

118.57 

Y 

1 

0.65 

6.23933 

0.026033 

0.802 

168.16 

30.05 

284.95 


2 

0.0737 

9.97199 

0.04723 

5.60 

348.18 

298.78 

49.35 

Z 

1 

0.168 

5.45138 

0.13299 

1.453 

133.88 

186.20 

113.99 


2 

0.934 

8.78422 

0.01305 

0.204 

313.83 

137.53 

146.93 

were undetected 

HD 128311 c has 



6. DISCUSSION 


been detected astrometrically (McArthur et al. 2014), 
but its companion planet has not, so the mutual Inciina- 
tion remains unknown. 


5.1. 2:1 Systems: HD 128311 and HD 73426 

We performed 100 integrations of HD 128311 and 
HD 73526 for 10 Myr each. We found no stable con¬ 
figurations of the former if the orbits were allowed to be 
non-planar. This should not be taken to mean the sys¬ 
tem must be copla nar as we may n ot h ave considered 
enough cases. [Rein & Papaloizou (2009) found the sys¬ 
tem could form in a co planar con figuration through con¬ 
vergent migration, and McArthur et al. (2014) found the 
best-fit coplanar solution to the system is dynamically 
stable and not in resonance. At this point, we conclude 
that this system is likely to be in a coplanar configu¬ 
ration, which precludes the chaotic evolution we report 
here. 

HD 73526, on the other hand, could be evolving chaot¬ 
ically. In Table [9] we list three versions that are stable for 
10 Gyr and show chaotic evolution. Note that we always 
use a stellar mass of 1.08 M e . In Fig. fl6| we show the 
evolution of System AB. Unlike the previous systems, the 
resonant arguments do not appear to switch between li- 
bration and circulation, and the i-arguments do not show 
any libration at all. The Q\ argument librates about 0 
for the full 10 Gyr year duration. Nonetheless, the ec¬ 
centricities and inclinations appear to be coupled and to 
switch between modes, as was seen in § [3] Systems AA 
and AC show similar behavior with similar amplitudes. 
Of the previous systems, HD 73526 is most similar to 
System B (c.f. Fig. [ 5 }. 


5.2. The 3:1 System HD 60532 

In this section we consider the HD 60532 system which 
is in the 3:1 MMR. In Table flOl we list 3 cases which show 
chaos for 10 Gyr, and present the evolution of case BC in 
Fig. 17 As with HD 75326, the i-resonance arguments 
circulate, but the e-arguments switch between libration 
and circulation. This system appears qualitatively sim¬ 
ilar to those in §§2- Over 10 Gyr, the system evolves 
chaotically, as shown in the right panels. The BA and BB 
systems are qualitatively similar, but the mode-switching 
is not as dramatic. 

5.3. The 3:2 System HD 45364 

We found no configurations of HD 45364 that were sta¬ 
ble for 10 Gyr, but one trial did survive for 4.557 Gyr 
while conserving ene rgy to 1 part in 10 6 . Its initial con¬ 
ditions are in Table 
Fig-fTSl 


The previous three section have demonstrated that 1) 
planets in an MMR and with mutual inclinations can 
experience chaotic evolution of orbital elements for Gyr 
while at least 1 resonant argument librates throughout, 
2) planet-planet scattering, in which one planet is re¬ 
moved from a planetary system by gravitational inter¬ 
actions, can leave behind two planets in the 2:1 MMR 
and with significant mutual inclinations, and 3) several 
systems known to be in an MMR could have mutual incli¬ 
nations that induce chaotic evolution, but maintain the 
resonance. In this section, we discuss the theoretical and 
observational implications of these results, as well as de¬ 
scribe the limits of our analysis, which naturally leads to 
directions for future research on this topic. 

Fig- HI shows that conjunction can librate about mul¬ 
tiple centers, some of which, such as vj\ , are predicted 
by classic celestial mechanics. These libration centers 
are derived from low-order e xpansions of the “ disturbing 
function” (for a review see Murray & Dermott (1999)). If 
one term dominates, there are specific stable longitudes 
for conjunction, e.g. w\ if the term 2 A 2 — Ai — vo\ domi¬ 
nates. However, with large and varying values of e and i 
multiple terms are important. Stable longitudes can mi¬ 
grate or become unstable, allowing transitions of libra¬ 
tion to different kinematic modes. This movement could 
be due to the deepening of nearby minima as e and/or 
i changes. As the orbits evolve, conjunction could gain 
access to another minimum, leading to a change in the 
libration center. 

The behavior has multiple drivers that each behave like 
a pendulum. In effect, the system is analogous to a com¬ 
pound pendulum, e.g. Eq. '©• As the system evolves, 
the planets are able to move into different modes of os¬ 
cillation. For some modes, the stable longitude is repre¬ 
sented by classic longitudes like w \, but others may only 
be derivable using higher order theory. Our hypothesis 
should be testable from derivation of high-order models 
of resonant behavior from the disturbing function. Such 
an analysis was beyond the scope of this work, which is 
just a demonstration of the amplitude and duration of 
the chaos, but is clearly desirable. 

This result has important implications for theoretical 
work on orbital stability. Common approaches for iden¬ 
tifying unstable orbits rely on short integrations (~ 1000 
orbits) and a subsequent analysis of the orbital evolution 
to count the number of frequencies in the orbital oscil¬ 
lations: a larger number could indicate the system is 
chaotic and unstable. Examples include the Mean Expo¬ 


ll and its evolution is shown in 


cotta & Simo 

2000 

' j - ,1 - - - _-T| ? Id 

Gozdziewski et ai. 2001)) and fre 

quency maps [e.g. 
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TABLE 9 

Initial Conditions for Selected HD 73526 Cases 


System 

Body 

m (M Jup ) 

a (AU) 

e 

H°) f 

nn f 

" (°) 

M(°) 

AA 

b 

2.89 

0.66857 

0.1728 

1.37 

182.66 

193.96 

96.27 


c 

2.576 

1.0615 

0.1126 

1.21 

2.44 

147.97 

107.05 

AB 

b 

2.881 

0.64347 

0.1824 

1.27 

262.07 

57.63 

95.97 


c 

2.405 

1.0311 

0.05134 

1.18 

81.74 

8.87 

92.51 

AC 

b 

2.87 

0.6572 

0.1552 

3.37 

342.84 

349.04 

75.29 


c 

2.691 

1.0274 

0.08992 

2.85 

162.79 

324.05 

77.07 


f Measured relative to invariable plane 


TABLE 10 

Initial Conditions for Selected HD 60532 Cases 

System 

Body 

m (M-Jup) 

a (AU) 

e 

i («>)/ 

« (°) / 

(°) 

M (°) 

BA 


b 0.997 

0.759 

0.3037 

2.56 

267.21 

185.73 

162.6 


c 

2.42 

1.5943 

0.1578 

0.692 

87.1 

252.05 

322.57 

BB 

b 

1.053 

0.7582 

0.2978 

1.403 

253.26 

212.28 

277.17 


c 

2.464 

1.5661 

0.0354 

0.399 

73.39 

250.17 

322.34 

BC 

b 

1.062 

0.7587 

0.2816 

5.806 

293.03 

79.61 

209.64 


c 

2.545 

1.572 

0.0207 

1.62 

113.05 

341.14 

321.88 


f Measured relative to invariable plane 


TABLE 11 

Initial Conditions for the HD 45364 Case 


Body 

m (Mjup) 

a (AU) 

e 

i(°) f 


(°) 

M(°) 

b 

0.1874 

0.6822 

0.1782 

2.07 

359.31 

110.52 

199.03 

c 

0.6581 

0.8969 

0.09935 

0.509 

179.38 

320.73 

265.9 


Measured relative to invariable plane 


Our results suggest that those methods are susceptible 
to labeling long-lived systems as short-lived. The inves¬ 
tigation of viable orbital architectures of planets in an 
inclined MMR should seek configurations that can per¬ 
sist for the age of the host star as described here, which 
appear unstable on the short term, but are in fact long- 
lived. For exa mple, HD 202206 is likely in a 5:1 MMR, 


but Correia et al. (2005) used a frequency analysis to con¬ 
strain the orbital parameters rather than N-body mod¬ 
els. Future work should explore the the veracity of these 
approximate methods for the case of inclined MMRs. 

At epochs of very high e, we expect some planets to 
tidally circularize. For example, Systems D and U could 
not persist as shown because tides would circularize the 
planet. Two scenarios are plausible, depending on the 
dissipation rate in the planet. If the dissipation is rela¬ 
tively weak, the resonance may be maintained such that 
the resonant pair migrates inward together. Such sys¬ 
tems should be represented in radial velocity and transit 
surveys, which are biased toward the detection of plan¬ 
ets on relatively short orbital periods. The three cases 
presented in § 5 are possible examples. In principle, the 
inward migration could be arrested if the tidal dissipa¬ 
tion forces the pair into a configuration that diminishes 
the maximum eccentricity. Thus, planets found far from 
the host star today could nonetheless have formed at 
larger distances and moved in. As the dissipation can be 
episodic, it could take a long time for the pair to migrate 
inward significantly, further increasing the likelihood to 


detect the planets where tidal forces are expected to be 
insignificant. Thi s possibility has been discussed for non - 
resonant systems (Li et al. 2014 Dawson & Chiang 2014), 
and we find it can also occur for MMRs. Future work 
should explore the role of tidal dissipation in these sys¬ 
tems and determine if there is any predicted or observed 
signature of tidal evolution resulting from an inclined 
MMR. 

On the other hand, if the tidal dissipation in the planet 
is very large, it may break the resonance, resulting in 
rapid tidal circularization and migration, leaving one 
close-in planet and one more distant a planet. Although 
many close-in planets are singletons (Steffen et al. |2012 ), 
some are known to host more distant companions . These 
companions have been suggested as perturbers that could 
drive eccent ricities to high values through Kozai-type 
oscillations (Takeda & Rasio 2005), but the possibility 
of resonant excitation of eccentricity has not previously 
been proposed. 

Note that these periods of high-e can occur when i 
has a wide range of values. Thus, the planet’s final or¬ 
bital plane could be independent of its initial orbital 
plane and be misaligned with its host st ar’s s pin axis. 
Such spin-orbit misalignm ent is observed (Triaud et al. 


2010 Hirano et al. 20141, and numerous mechanisms 


have been proposed, such as pla net-planet scattering fol - 


lowed by tidal circularization (Chatterjee et al. 2008) 
tidal circularization during phases of high eccent ricity 
due to interactions with very distant perturbers (Fab- 
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| Fig. 14.— Orbit al evolution of System W in the same format as Fig. [8] 
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| Fig. 15.— Orbital evolution of System X in the same format as Fig. [8] 
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| Fig. 16.— Orbital evolution of System AB (HD 73526) in the same format as Fig, [s] 
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rycky fe Tremaine|2007 Storch et al.|2014 |, and interac- 
tions wi th the gas disk and a distant stellar binary com¬ 
panion (Batygin| 2012). The chaotic evolution of e and i 
in an inclined MMR is another process to add to this list. 
However, given the low occurrence rate of chaotic incli¬ 
nation resonances formed by scattering, it is unlikely to 
be a dominant mechanism. But note that this inference 
is based on only one set of scattering simulations with 
specific initial conditions - others may be more efficient 
and producing these types of systems. 

Inclined MMRs could impact the distribution of period 
ratios of planets found via transit by making one planet’s 
orbital plane significantly different from the other. After 
applyi ng a geometric transit correction, [Steffen & Hwang] 
(2014) find a relative excess of planets in the 3:2 reso- 


The most likely instru ment to discover planet s in an in- 


2008 


8 J), as shown 
d in § [5] are 


nance and relative deficits in 2:1 and 3:1 in Kepler data. 
However, for the systems in which only 2 planets are de¬ 
tected, the excess of 3:2 pairs disappears, the 2:1 is still 
depleted, and the 3:1 appears to be unaffected. If addi¬ 
tional planets destabilize inclined MMRs, then the trend 
in the 3:2 MMR may be indicative of inclined MMRs, 
as we would not expect systems that evolve with large 
amplitude to have companions nearby. However, caution 
is necessary when interpreting these data, as knowledge 
of the underlying population and the role of tidal evo¬ 
lution is critical. Migration during the protopla netary 


disk phase often leads to capture into resonance (Snell- 


grove et al. 2001 Lee & Peale 
mordial excess of planets in IVilv 


2002 ), producing a pri- 
Rs. On the other hand, 


tidal evolution of planets in MMRs tends to pull them to 


period ratios that are not exactly commensurate (Lith- 


clined MMR is GAIA (Casertano et al. 
in Fig. [9] The known systems we examine 
all goodicandidates for GAIA astrometry, and so in the 
next 5 years we should find out if any of them could be 
experiencing chaotic evolution. 

We note that our simulations of known systems failed 
to reveal any in which the i arguments librate. Instead 
the e-arguments switch between circulation and libra- 
tion, which likely drives the chaos. Note that System B 
evolves in a similar manner, and many other cases do 
as well. Of our 3 systems formed b y sc attering, 2 were 
similar to the known systems (Fig. [l5| ). A further ex¬ 
ploration of known systems, either by including more or 
broadening the parameter space survey, could reveal if 
this trend is real. If known systems are only consistent 
with i-argument circulation, it could provide important 
clues to the formation and frequency of exoplanets in 
inclined MMRs. 

Throughout this study we have described systems that 
survive for 10 Gyr as “stable.” As most of our host stars 
are solar analogs, this usage is reasonable as the systems 
survive for the main sequence lifetimes of the host stars. 
However, in some cases we find destabilization can occur 
on > 1 Gyr timescales. Hence, these systems are not 
“stable” in the sense that they could survive indefinitely. 
Late-term destabilization of systems in inclined MMRs 
may explain the observation that the host stars of planets 
in the 2:1 MMR tend to be younger than other host stars 
(Koriski & Zucker[2011). If this observation is validated, 
it would support the hypothesis that some MMRs evolve 


wick & Wu 

2012; Batygin & Morbidelli 2013] 

Delisle & 

chaotically and disrupt on long timescales. 

Laskar 201- 

). Moreover, the formation of inclined reso- 

In the previous sections we did not consider the role 


be intrinsically rare. Thus, it is not obvious that inclined 
MMRs are sculpting the close-in exoplanet population, 
but our results suggest that they could. 

Future work should identify boundaries to the long- 
lived but chaotic resonances discovered here and employ 
a quantitative description of the chaos. All of our sim¬ 
ulations of hypothetical systems began with planetary 
orbital periods at exact commensurability, but this state 
is not typically observed for exoplanets, see Tableland 
§ 5. A large suite of N-body simulations that integrate 
systems to 10 Gyr is required to map out the bound¬ 
aries of the chaotic and resonant behavior. Throughout 
this study we have referred to systems as chaotic as they 
clearly are. However, as a system moves toward regular 
motion, the motion might no longer be obviously chaotic. 
The use of Lyapunov exponents or other metrics would 
be necessary to elucidate the true boundaries of the long- 
lived chaotic motion. 

Identifying inclined MMRs in transit and /or RV data is 
possible but difficult (see e.g. Dawson et al.||2014 ), but 
astrometric measurements are probably the best route 
to find their existence. HST has successfully detected 
some plane ts for which RV detections suggest an MMR 


longevity of a chaotically evolving system. It is clear that 
for system in which e\ ~ 1 that interior planets are for¬ 
bidden. However, exterior planets could exist provided 
they are distant and/or small. On the other hand, some 
systems show low amplitude chaotic variations and may 
therefore be robust to perturbations from a third planet. 
We do not map out the role of a third companion in the 
stability of our systems, but future work should explore 
how they modify the results presented here. 

Most of our simulations begin with an Earth-mass 
planet at 1 AU from a solar-m ass star, and would be con- 


parapu et al. 


2013). If habitable, these worlds would be 


sidered potentially habitable (|Kasting et al. 1993 Kop- 

TIs 


markedly different from the Earth, with unpredictable 
climates on geologic timescales. For planets in which 
the eccentricity grows large, the planet could occasion¬ 
ally enter a runaway greenhouse (incident radiation flux 
scales as (1—e 2 ) -1 / 2 ) rendering the planet uninhabitable. 
Large inclination fluctuations can also lead to large vari¬ 
ations in obliquity, which is a major driv er of climate 
evolu tion (e.g. Williams fe Kastiiig||1997 Spiegel et al. 


is present ( 

Benedict et al. 

2002 

McArthur et al. 20141 

and two planets not in an MMR 

McArthur et ai.|2U10), 


planet of HD 128311 has been detected, and the in¬ 


ner is potentially detectable with HST (|McArthur et al. 
20141, but its precarious architecture relative to dynam- 


ical instability suggests it is in a coplanar configuration. 


2009). While extremely fast and large variations of obliq- 
uity could be detrimental to habitability, at the outer 
edge of the habitable zone, these variations can suppress 
ice sheet gro wth and, in principle, increase a planet’s 
habitability (Armstrong et al. 2014). For planets that 
occasionally reach very high eccentricity, tidal dissipa¬ 
tion could ult imatel y pull the planet out of the habitable 
zone (Barnes et al.][20081. Should any potentially habit- 
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Fig. 18. — Orbital evolution of the HD 45364 case in the same format 
able planets be found in an MMR, it will be imperative 
to understand the orbital evolution, and its connection to 
climate, prior to investing spectroscopic observati ons of 


its atmosphere to search for biosignatures (see e.g. Peril¬ 
ing et al. 1120091 |Kaltenegger & Traub||2009[ |lVlisra et al. 
2014T 


7. CONCLUSIONS 

We have simulated the orbital evolution of exoplanets 
in mean motion resonances and inclinations and found 
the orbits can evolve chaotically for at least 10 Gyr. We 
hypothesize that these systems behave like compound 
pendula, which are naturally chaotic systems that can 
switch between modes of oscillation, as seen in our sim¬ 
ulations (see Fig. [3]). We find this chaotic motion over a 
range of mass ratios and for the 2:1, 3:2 and 3:1 res¬ 
onance. We also tested different N-body codes using 
different integration schemes, and conclude the results 
are robust. Inclined MMRs can be produced by planet- 
planet scattering and the resultant systems are qualita¬ 
tively similar to our simulations of known systems in an 
MMR. 

These results have numerous implications for both the¬ 
ory and observations: 

• Approximate methods to estimate stability with 
short integrations may be unreliable near MMRs. 

• Close-in planets may arrive at their current orbits 


as Fig. Is] 

clue to eccentricity excitation by inclined MMRs 
followed by rapid tidal circularization. 

• Some short-period planets with orbital planes mis¬ 
aligned with the stellar spin axis may be produced 
by systems initially in inclined MMRs. 

• MMR pairs may be episodically migrating inward 
due to weak dissipation occurring during epochs of 
very large eccentricity. 

• Systems in an MMR may be systematically 
younger than other multiplanet systems due to the 
destabilization of older MMR systems. 

• The distribution of period ratios of adjacent plan¬ 
ets detected via transit may be skewed by inclined 
MMRs. 

• Potentially habitable planets may be severely im¬ 
pacted by the orbital architecture of the system. 

Although no systems are currently known to demon¬ 
strate the behavior we have outlined here, the GAIA 
space telescope has the power to detect hundreds of giant 
exoplanets in inclined MMRs. 

This work was supported by NASA’s Virtual Plan¬ 
etary Laboratory under Cooperative Agreement No. 
NNA13AA93A and NSF grant AST-1108882. 
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